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Abstract 



We study the growth of correlations in systems with weak long-range interactions. 
. Starting from the BBGKY hierarchy, we determine the evolution of the two-body corre- 

Ch I lation function by using an expansion of the solutions of the hierarchy in powers of 

in a proper thermodynamic limit — > -|-oo. These correlations are responsible for the 
"collisional" evolution of the system beyond the Vlasov regime due to finite N effects. We 
I obtain a general kinetic equation that can be applied to spatially inhomogeneous systems 

K*" ' and that takes into account memory effects. These peculiarities are specific to systems 

■ with unshielded long-range interactions. For spatially homogeneous systems with short 

, mcmory time like plasmas, we recover the classical Landau (or Lenard-Balescu) equations, 

'sj" I An interest of our approach is to develop a formalism that remains in physical space (in- 

' stead of Fourier space) and that can deal with spatially inhomogeneous systems. This 

, enlightens the basic physics and provides novel kinetic equations with a clear physical 

I interpretation. However, unless we restrict ourselves to spatially homogeneous systems, 

closed kinetic equations can be obtained only if we ignore some collective effects between 
particles. General exact coupled equations taking into account collective effects are also 
given. We use this kinetic theory to discuss the processes of violent collisionless relaxation 
' and slow collisional relaxation in systems with weak long-range interactions. In particular, 

we investigate the dependence of the relaxation time with the system size and provide a 
coherent discussion of all the numerical results obtained for these systems. 



1 Introduction 

Systems with long-range interactions are numerous in nature pLj. Some examples include self- 
gravitating systems, two-dimensional vortices, neutral and non-neutral plasmas, bacterial pop- 
ulations, defects in solids, etc... When the potential of interaction is attractive and unshielded, 
these systems can spontaneously organize into coherent structures accounting for the diversity 
of the objects observed in the universe. For example, self-gravitating systems organize into 
planets, stars, galaxies, clusters of galaxies... On the other hand, two-dimensional turbulent 
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flows organize into jets (like the gulf stream on the earth) or large-scale vortices (like Jupiter's 
great red spot in the jovian atmosphere). Biological populations (like bacteria, amoebae, en- 
dothelial cells,...) also interact via long-range signals through the phenomenon of chemotaxis. 
Chemotactic aggregation leads to the spontaneous appearance of patterns hke stripes and spots, 
filaments, vasculature,... Although these astrophysical, hydrodynamical and biological systems 
are physically different, they share a lot of analogies due to the long-range attractive nature of 
the potential of interaction [2]. 

In view of the complexity of these systems, it is natural to try to understand their struc- 
ture and organization in terms of statistical mechanics pj. Since systems with long-range 
interactions are generically spatially inhomogeneous, it is clear at first sights that the usual 
thermodynamic limit +00 with N/V fixed is not valid. Therefore, the ordinary methods 
of statistical mechanics and kinetic theory must be reformulated and adapted to these systems. 
We shall assume, however, that the basic concepts are not altered so that the description of 
these systems must be done in consistency with the foundations of statistical mechanics and 
kinetic theory. In previous papers of this series 0, S] (denoted Papers I and II), we have un- 
dertaken a systematic study of the dynamics and thermodynamics of systems with long-range 
interactions. In Paper I, we have considered the statistical equilibrium states and the static 
correlation functions. We have shown that there exists a critical temperature Tc (for Brownian 
systems) or a critical energy Ec (for Hamiltonian systems) above which the system is spatially 
homogeneous and below which the homogeneous phase becomes unstable and is replaced by a 
clustered phase. In Paper II, using an analogy with plasma physics, we have developed a kinetic 
theory of systems with long-range interactions in the homogeneous phase. In the present paper 
(Paper III), we propose new derivations of the kinetic equations that take into account non- 
markovian effects and that can be applied to spatially inhomogeneous configurations. These 
extensions are specific to systems with unshielded long-range interactions and they are novel 
with respect to the much more studied case of neutral plasmas. They complete the results of 
Paper II that were only valid for spatially homogeneous and markovian systems. However, a 
limitation of the present approach is to neglect collective effects. These effects were taken into 
account in Paper II for spatially homogeneous systems. 

This paper is organized as follows. In Sec. [21 we derive a general kinetic equation for 
Hamiltonian systems with weak long-range interactions from the BBGKY hierarchy. This 
equation is valid at order 0{1/N) in an expansion of the solutions of the equations of the 
hierarchy in powers of in the proper thermodynamic limit — > +00 defined in Paper I. 
For A^ — i> +00, this kinetic equation reduces to the Vlasov equation. At order 0{1/N) it takes 
into account the effect of "collisions" (more properly "correlations") between particles due to 
finite A^ effects (graininess) . It describes therefore the evolution of the system on a timescale 
A^t^i, where ti) is the dynamical time. This general kinetic equation applies to systems that 
can be spatially inhomogeneous and takes into account non-markovian effects. If we restrict 
ourselves to spatially homogeneous systems and neglect memory terms, we recover the Landau 
equation as a special case. In Sees. [3] and HI we use this kinetic theory to discuss the processes 
of violent collisionless relaxation and slow coUisional relaxation in systems with weak long-range 
interactions. We review several results obtained for self-gravitating systems, two-dimensional 
vortices and the HMF model, emphasize their connections and try to explain them in the 
light of the kinetic theory. In particular, we investigate the dependence of the relaxation time 
with the system size. We also propose a scenario according to which, for a large class of 
initial conditions, the transient states of the coUisional relaxation of the HMF model could be 
described by spatially homogeneous Tsallis distributions (polytropes) with a compact support 
and with an index q{t) > 1 slowly decreasing with time until they become Vlasov unstable and 
relax towards the Boltzmann distribution. 
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2 Kinetic equation from the BBGKY hierarchy 



In this section, we derive a general kinetic equation fl33|) for Hamiltonian systems with weak 
long-range interactions. We start from the BBGKY hierarchy and use a systematic expansion 
of the solutions of the equations of this hierarchy in powers of in a proper thermodynamic 
limit N ^ +00. The kinetic equation (l33l) is valid at order 0{1/N). 



2.1 The expansion 

We consider a system of particles with long-range interactions described by the Hamiltonian 
equations (I-l). Basically, the evolution of the A^-body distribution function is governed by 
the Liouville equation (1-2). Introducing the reduced probability distributions (1-6), we can 
construct the complete BBGKY hierarchy (II-l). The first two equations of this hierarchy, 
governing the evolution of the one and two-body distributions Pi(xi,t) and P2(xi,X2,t), are 
given by Eqs. (II-2) and (II-3). We recall that x stands for (r, v). We now decompose the 
distributions functions in the form (1-14) and (1-15) where P2(xi,X2,t) and P3(xi, X2, X3, t) are 
the two and three-body correlation functions (or cumulants). Substituting this decomposition 
in Eq. (II-2), we first obtain 

■ + (AT - 1)^ / F(2 ^ l)Pi(x2)c/x2 



dt dri dvi 

+ (Ar-l)Ay F(2^1)P^(xi,X2)dx2 = 0, (1) 

where F(j i) is the force by unit of mass created by particle j on particle i. It is related to 
the potential of interaction Uij = u{\ri — rj\) by 

F(J -> = -J-^. (2) 

Then, substituting the decompositions (1-14) and (1-15) in (II-3) and using Eq. ([1]) to simplify 
some terms, we get 

-P,i^x,)—j F(3 l)Pi(xi)Pi(x3)rfx3 

d 



F(3^1)P^(xi,X3)Pi(x2)dX3 



+ {N-2)—J F(3^1)P^(Xi,X2)Pi(x3)rfX3 
+ {N-2)—J F(3^1)P^(x2,X3)Pl(xi)rfX3 

+ {N - 2)— / F(3 ^ l)P3(xi, X2, X3)rfx3 + (1 ^ 2) = 0. (3) 

Equations ([I])-([3]) are exact for all N but the hierarchy is not closed. We shall now consider 
the thermodynamic limit defined in Paper I. It corresponds to ^ +00 in such a way that 



^In Paper II, some terms were missing in Eq. (II-5) of the BBGKY hierarchy because we systematically took 
— 1 ~ iV and — 2 ~ which is not correct if we consider terms of order 0{1/N). 
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the normalized temperature r] = jSNrn^u^: and the normahzed energy e = E /{u^N'^im?) are 
fixed, where represents the typical strength of the potential of interaction. In general, 
the potential of interaction is written as u{vij) = ku{vij) where k is the coupling constant 
(e.g., G for self-gravitating systems or k for the HMF model). By a suitable normalization 
of the parameters, this thermodynamic limit is such that the coupling constant behaves like 
k ^ Ut, ^ 1/N while the individual mass m ~ 1, the inverse temperature /3 ~ 1, the energy 
per particle E/N ~ 1 and the volume V ^ 1 are of order unity q This implies that |x| ~ 1 
and |F(j — >• i)\ ~ On the other hand, the dynamical time tr, ~ R/vtyp ~ l/y/kp ~ 1 

is of order unity (p ~ M/V is the average density and the typical velocity Vtyp has been 
obtained by equating the kinetic energy ~ Nmv'^ and the potential energy ~ N'^m^ku). Since 
the normalized coupling constant f3m'^u^ = rj/N ~ goes to zero for — * +oo, we are 
studying systems with weak long-range interactions. It is argued in Papers I and II that there 
exists solutions of the whole BBGKY hierarchy such that the correlation functions scale 
like 1/N^~^. This implicitly assumes that the initial condition has no correlation, or that the 
initial correlations respect this scaling (if there are strong correlations in the initial state, the 
system will take a long time to erase them and the kinetic theory will be different from the 
one developed in the sequel). If this scaling is satisfied, we can consider an expansion of the 
solutions of the equations of the hierarchy in terms of the small parameter This is similar 
to the expansion in terms of the plasma parameter made in plasma physics. However, in plasma 
physics the systems are spatially homogeneous while, in the present case, we shall take into 
account spatial inhomogeneity. This brings additional terms in the kinetic equations that are 
absent in plasma physics. Therefore, strictly speaking, the hierarchy that we consider is different 
from the ordinary BBGKY hierarchy. Recalling that Pi ~ 1, ~ l^U — ^ 0I ~ 1/^, 

we obtain at order 



dPi 
~dt 



^i— h (A/ - 1) 



dri 



+N 



d 



F(2 - 



F(2^ l)Pi(x2)cix2 

l)P^(xi,X2)rfX2 = 0, 



(4) 



dPL dPi 
at or I 



F(2 ^ 1) - / F(3 ^ l)Pi(x3)rfx3 



^l(x2)^(Xi, 



dP' C d C 

-fiV^ / F(3 ^ l)Pi(x3)rfx3 + N— I F(3 1)P^(X2, X3)Pi(xi)rfx3 + (1 ^ 2) = 0. (5) 

If we introduce the notations / = NmPi (distribution function) and g = iV^Pg (two-body 
correlation function), we get 



dt OTi N ovi 



-m- 



d_ 



F(2^1)(7(xi,X2)dx2, 



(6) 



d 



| + v,|»+(F),|^ + V(2^1)a|A 

ot OTi avi avi 

- /F(3^1)(7(x2,X3,t)-rfx3 + (1^2) = 0, 
1 J ^ 



(7) 



^Alternatively, we can assume that the mass of the particles scales like m ~ while /c ~ ~ 1, /3 ~ TV, 
E ^ 1 and V ^ I. In this scaling, the total mass M Nra is of order unity. 



4 



where we have introduced the abbreviations /i = /(ri,vi,t) and /2 = /(r2,V2,t). We have 
also introduced the mean force (by unit of mass) created in ri by all the particles 



(F)i = /f(2 ^ I)^dr2rfv2 
J rn 



and the fluctuating force (by unit of mass) created by particle 2 on particle 1: 

^(2^1) = F(2^1)-l(F)i. (9) 

These equations are exact at the order 0{1/N). They form therefore the right basis to develop 
a kinetic theory for Hamiltonian systems with weak long-range interactions. We note that 
these equations are similar to the BBGKY hierarchy of plasma physics but not identical. One 
difference is the (A^ — 1)/N term in Eq. (Q. The other difference is the presence of the 
fluctuating force JF(2 1) instead of F{2 1) due to the spatial inhomogeneity of the 
system. In plasma physics, the system is homogeneous over distances of the order of the Debye 
length so the mean force (F) vanishes. 

2.2 The Vlasov equation and beyond 

Recalling that P2 ^ 1/N, we note that 

P2(X1, X2, t) = Pi(xi, t)Pi(x2, t) + 0{1/N). (10) 

If we consider the limit — > +00 (for a fixed time t), we see that the correlations between 
particles can be neglected so that the two-body distribution function factorizes in two one- 
body distribution functions i.e. P2(xi,X2,t) = Pi(xi, t)Pi(x2, t). Therefore the mean field 
approximation is exact in the limit —>■ +00. Substituting this result in Eq. ([1]), we obtain 
the Vlasov equation 

f + v,|i + (F).|i = 0. (11) 

at OTi avi 

This equation describes the collisionless evolution of the system up to a time at least of order 
Ntj:) (where to is the dynamical time). In practice, iV S> 1 so that the domain of validity of 
the Vlasov equation is huge (for example, in typical stellar systems ~ 10^ — 10^^). When the 
Vlasov equation is coupled to an attractive unshielded long-range potential of interaction, it 
can develop a process of violent relaxation towards a quasi stationary state (QSS). This process 
will be discussed specifically in Sec. El 

If we want to describe the collisional evolution of the system, we need to consider finite A^ 
effects. Equations ([6]) and ([7]) describe the evolution of the system on a timescale of order Nto- 
The equation for the evolution of the smooth distribution function is of the form 

+ "^^^^ + -ir^^^'d;r^ - ^^t^^' ^^^^ 

where C^v is a "collision" term analogous to the one arising in the Boltzmann equation. In the 
present context, there are not real collisions between particles. The term on the right hand side 
of Eq. (fT2!) is due to the development of correlations between particles as time goes on. It is 
related to the two-body correlation function (7(xi, X2, t) which is itself related to the distribution 
function /(xi, t) by Eq. Our aim is to obtain an expression for the collision term CN[f] at 
the order 1/A^. The difficulty with Eq. (JTj) for the two-body correlation function is that it is an 
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integrodifferential equation. The second term is an advective term, the third term is the source 
of the correlation and the third term takes into account the retroaction of the system as a whole 
due to a change of the correlation function. In this paper, we shall neglect the contribution of 
the integral in Eq. ([7]). Then, we get the coupled system 



dt dri 



+ 



N 



N 



-m- 



d_ 



F(2^1)(7(xi,X2)rfx2, 



(13) 



do 

— + 

dt 



d d d 

Vl7^+V2— + F 1 — 

OTi OT2 avi 



_d_ 

dV2 



+ 



T{2 



I — 
9vi 



Til 



I — 

(9V2 



m m 



(14) 



The integral that we have neglected contains "collective effects" that describe the polarization of 
the medium. In plasma physics, they are responsible for the Debye shielding, i.e. the fact that a 
charge is surrounded by a polarization cloud of opposite charges that diminish the interaction. 
These collective effects are taken into account in the Lenard-Balescu equation through the 
dielectric function (see Paper II) . However, this equation is restricted to spatially homogeneous 
systems and based on a Markovian approximation. These assumptions are necessary to use 
Laplace- Fourier transforms in order to solve the integro-differential equation (jT}). Here, we 
want to describe more general situations where the interaction is not shielded so that the system 
can be spatially inhomogeneous. If we neglect collective effects, we can obtain a general kinetic 
equation in a closed form that is valid for systems that are not necessarily homogeneous and 
that can take into account memory effects. This equation has interest in its own right (despite 
its limitations) because its structure bears a lot of physical significance. Before deriving this 
general equation, we shall first consider the case of spatially homogeneous systems and make 
the link with the familiar Landau equation. 



2.3 The Landau equation 

For a spatially homogeneous system, the distribution function and the two-body correlation 
function can be written / = /(vi,t) and g = (^(vi, V2,ri — r2,t). In that case, Eqs. f|T3l) -f|T^ 
become 



dfi ^ d f du 



(15) 



dg dg du f d 

h w ■ — = — • 

dt (9x (9x V9vi dv2 



^ \(vi,t)^(v2,t), 

m 



(16) 



where we have used the fact that F(l ^ 2) = — F(2 1) and noted x = ri— r2 and w = vi — V2. 
Taking the Fourier transform of Eq. ( |T6l) and introducing the notations d = d/d\i — d/dv2-, 
fi = /(vi,t) and /2 = /(v2,t), we obtain 



+ ik-wg = —u{k)k ■ S/i/a. 



dt m 



(17) 



In terms of the Fourier transform of the correlation function, the kinetic equation f[T^ can be 
rewritten 



k'u(fc)Im^(vi, V2, k, t) (ik(iv2. 
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We shall assume that Im^(vi, V2, k, t) relaxes on a timescale that is much smaller than the 
timescale on which f{vi,t) changes. This is the equivalent of the Bogoliubov hypothesis in 
plasma physics. If we ignore memory effects, we can integrate the first order differential equation 
(fT7|l by considering the last term as a constant. This yields 



^(vi,V2,k,t)= f dT-ku{k)e-''^-'"^dh{t)f2{t), (19) 
Jo rn 

where we have assumed that no correlation is present initially: g{t = 0) = 0. Then, we can 
replace Im5f(vi, V2, k, t) in Eq. (ITS!) by its value obtained for t +oo, which reads 

Im^(k,vi,V2,+oo) = -ku{k)6{k-w)dfi{t)f2{t). (20) 

m 

Substituting this relation in Eq. ( ITSj) . we obtain the Landau equation in the form 

f = ^P')*"s!f / "^^"^'^^'^-m^H^ ■ w) (/.|| -/.§). (21) 

Other equivalent expressions of the Landau equation are given in Paper II. The Landau equation 
ignores collective effects. Collective effects can be taken into account by keeping the contri- 
bution of the last integral in Eq. (^^. For spatially homogeneous systems, the calculations 
can be carried out explicitly in the complex plane [5J and lead to the Lenard-Balescu equation 
discussed in Paper II (the Lenard-Balescu equation can be obtained from the Landau equa- 
tion by replacing the potential u{k) by the "screened" potential u{k)/\e(k,'k ■ V2)| including 
the dielectric function). The Landau and Lenard-Balescu equations conserve mass and energy 
(reducing to the kinetic energy for a spatially homogeneous system) and monotonically increase 
the Boltzmann entropy f6] . The collisional evolution is due to a condition of resonance between 
the particles orbits. For homogeneous systems, the condition of resonance encapsulated in the 
5-function appearing in the Landau and Lenard-Balescu equations corresponds to k • vi = k ■ V2 
with vi 7^ V2. For d > 1, the only stationary solution is the Maxwell distribution. Because of 
the if-theorem, the Landau and Lenard-Balescu equations relax towards the Maxwell distribu- 
tion. Since the collision term in Eq. (^T^ is valid at order 0{1/N), the relaxation time scales 
hke 

tn ~ NtD, {d > 1) (22) 

as can be seen directly from Eq. fl2Tl) by dimensional analysis (comparing the l.h.s. and the 
r.h.s., we have l/t/j ~ ulN ~ 1/A^ while ~ R/vtyp ~ 1 with the scalings introduced in 
Sec. 12.11) . A more precise estimate of the relaxation time is given in [7]. For one-dimensional 
systems, like the HMF model, the situation is different. For d = 1, the kinetic equation fl2Tl) 
reduces to 



(23) 



Therefore, the collision term CnU'] vanishes at the order 1/A^ because there is no resonance. 
The kinetic equation reduces to df/dt = so that the distribution function does not evolve 
at all on a timescale ~ Nto- This implies that, for one-dimensional homogeneous systems, the 
relaxation time to statistical equilibrium is larger than Nt^. Thus, we expect that 

tn>NtD, id = l). (24) 

The fact that the Lenard-Balescu collision term vanishes in ID is known for a long time in 
plasma physics (see, e.g., the last paragraph in [8]) and has been rediscovered recently in the 
context of the HMF model P H [TU]. 
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2.4 The non Markovian kinetic equation 

The above kinetic equations rely on the assumption that the correlation function relaxes much 
more rapidly than the distribution function. The Markovian approximation is expected to be 
a good approximation in the limit — > +00 that we consider since the distribution function 
changes on a slow timescale of order NtD (where tu is the dynamical time) or even larger. 
However, for systems with long-range interactions, there are situations where the decorrelation 
time of the fluctuations can be very long so that the Markovian approximation may not be 
completely justified. This concerns in particular the case of self-gravitating systems for which 
the temporal correlation of the force decreases like 1/t (see [H] and Paper II). This is also the 
case for systems that are close to the critical point since the exponential relaxation time of the 
correlations diverges for E ^ 01 T ^ T^. (see [T2l [TO] and Paper II). Therefore, it can be of 
interest to derive non-markovian kinetic equations that may be relevant to such systems. If we 
keep the time variation of /(v,t) in Eq. (fT7|) . we obtain after integration 

^(vi,V2,k,t) = f dT-ku{k)e-''^'^^dMt-r)Mt-T). (25) 
Jo ^ 

Inserting this relation in Eq. f[T^ . we obtain a non Markovian kinetic equation 

^ = ^2^)'^^ / / d^^dkk^k'^uikf cos(k . wr) A _ l^j /(v„ t - r)f{^r„ t-r). 

(26) 

In particular, for the HMF model, using the notations of Paper I, we get 

We note that, when memory terms are taken into account, the collision term does not vanish. 
However, if we make the Markovian approximation /(fi,t — r) ~ f{vi,t), f{v2,t — T) ~ f{v2,t) 
and extend the time integral to infinity [^1, we obtain 

When memory terms are neglected we recover the fact that the Landau collision term vanishes 
for a spatially homogeneous one-dimensional system. By working close to the critical point in 
the HMF model (where the exponential relaxation time of the correlations diverges), it may 
be possible to see non-markovian effects in numerical simulations of the A^-body system. They 
should induce a small evolution of the homogeneous system on a timescale Nt^ as described 
by Eq. (l27|l or, more precisely, by its generalization taking into account collective effects (see 
Appendix R|) . This should not lead, however, to statistical equilibrium since Eq. (127|) clearly 
does not tend to the Boltzmann distribution. 



2.5 The kinetic equation for spatially inhomogeneous systems 



Relaxing the assumption that the system is spatially homogeneous, the equation (JT^ for the 
correlation function can be written 



dt 



T{2 



r 



_d_ 
(9vi 



T{1 



_d_ 



— (Xi,t) — (X2,t), 

m m 



(29) 



■^We could also consider an approximation where we make the Markovian approximation but keep the time 
integral going from to t (see Appendix . 
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where we have denoted the advective term by L (Liouvilhan operator) 
equation with the Green function 



= exp 



Solving formally this 



(30) 



we obtain 



C/(xi,X2,t) 



dTG{t,t- t) 



T{2 



d 



+ ^(1 



_d_ 



— (xi,t - r) — (x2,t - r). 
m m 



(31) 



The Green function constructed with the smooth field (F) means that, in order to evaluate 
the time integral in Eq. fl3T|) . we must move the coordinates rj(t — r) and Vj(t — r) of the 
particles with the mean field flow in phase space, adopting a Lagrangian point of view. Thus, 
in evaluating the time integral, the coordinates and Vj placed after the Greenian must be 
viewed as ri{t — r) and \i{t — r) where 



T 



Vi(t-s)rfs, Vi(t-r) = Vi(t) - / (F)(ri(t-s),t-s)rfs. (32) 



Substituting Eq. ([3U]) in Eq. ([I3D, we get 

dt dri N ^ ^^9vi 9< 



dr / dr2rfv2F^(2 l,t)G{t,t 



X 



d 



_d_ 



/(ri,vi,t 



/ 

r) — (r2,V2,t 
m 



T 



r). 



(33) 



This returns the general kinetic equation obtained by Kandrup [13] with the projection operator 
formalism (note that we can replace F'^[2 — > 1, t) by JF'^(2 — > l,t) in the first term of the r.h.s. 
of the equation since the fluctuations vanish in average). Equation fl5^ slightly differs from the 
equation obtained in |T3] by a term {N — 1)/N in the l.h.s. This new derivation of the kinetic 
equation (1331) from a systematic expansion of the solutions of the BBGKY hierarchy in powers 
of 1/A^ is valuable because the present formalism is considerably simpler than the projection 
operator formalism and clearly shows which terms have been neglected in the derivation. It 
also clearly shows that the kinetic equation l[3^) is valid at order 1/N so that it describes the 
"collisional" evolution of the system on a timescale of order Nto- 



2.6 Summary of the different kinetic equations 

Let us briefly summarize the different kinetic equations that appeared in our analysis. When 
collective terms are ignored, the kinetic equation describing the evolution of the system as a 
whole at order 1/A^ is 



df df N-1 
dt dr N ^ ' 



% = ^jy''l drid^r,F^{l 0)G(t, t - r) 



X {-^^(1 ^ + -^^(0 ^)^}^(^' V, t - r)^(r„ VI, t - r). (34) 

If we make a Markov approximation and extend the time integral to infinity, we obtain 
df df N-1 df d r^°° r 

i + + ^(F)^ ^ - drj dr,d.,F^il 0)G(t, t - r) 
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As we have indicated, the Markov approximation is justified for — > +00 so that Tcorr 
Ueiax ^ Nto- We do not assume, however, that the decorrelation time is "extremely" short 
(i.e, Tcorr ^ 0). Therefore, in the time integral, the distribution functions must be evaluated 
at (r(t — r), v(t — r)) and (ri(t — r), Vi(t — r)) where 

Vi{t-T) = r,{t)- Wi(t-s)ds, Vi(t-r) = Vi(t) - / {Y){ri{t - s),t)ds. (36) 

JQ Jo 

Comparing Eq. ( l36i) with Eq. ( 132|) . we have assumed that the mean force (F)(r,t) does not 
change substantially on the timescale Tcorr on which the time integral has essential contribution. 

For a spatially homogeneous system, using the fact that Eq. (132|) reduces to v(t — r) = v(t) 
since (F) = 0, Eq. takes the form 

If we make the integration on ri, using the relation (A3) of Paper II, we obtain the non 
markovian equation (|26l) . If we make the Markovian approximation /(v,t — r) ~ /(v,t), 
/(vi,t — r) ~ /(vi,t) and extend the time integration to +cx3, we get 

The integrals on r and ri can be performed as in Appendix A of Paper II and we finally obtain 
the Landau equation ( 12T]) . 

3 Violent collisionless relaxation 

In this section, we physically discuss the process of violent collisionless relaxation in relation 
with the Vlasov equation (fTTj) and point out several analogies between stellar systems, two- 
dimensional turbulence and the HMF model [I]. 

3.1 Quasi Stationary States 

When the Vlasov equation is coupled to an attractive unshielded long-range potential of inter- 
action it can develop a process of phase mixing and violent relaxation leading to the formation 
of a quasi-stationary state (QSS). This purely mean field process takes place on a very short 
time scale, of the order of a few dynamical times. This corresponds to the formation of galaxies 
in astrophysics, jets and vortices in geophysical and astrophysical fiows and clusters in the 
HMF model. Lynden-Bell [14j has proposed to describe these QSS in terms of statistical me- 
chanics, adapting the usual Boltzmann procedure so as to take into account the specificities 
of the Vlasov equation (in particular the conservation of the infinite class of Casimirs) |15]. 
This approach rests on the assumption that the collisionless mixing is efficient and that the 
ergodic hypothesis which sustains the statistical theory is fulfilled. There are situations where 
the Lynden-Bell prediction works relatively well. However, there are other situations where 
the Lynden-Bell prediction fails. It has been understood since the beginning [14J that violent 
relaxation may be incomplete in certain cases so that the Lynden-Bell mixing entropy is not 
maximized in the whole available phase space. Incomplete relaxation [16] can lead to more 
or less severe deviations from the Lynden-Bell statistics. Physically, the system tries to reach 
the Lynden-Bell maximum entropy state during violent relaxation but, in some cases, it can- 
not attain it because the variations of the potential, that are the engine of the evolution, die 
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away before the relaxation process is complete (there may be other reasons for incomplete 
relaxation). Since the Vlasov equation admits an infinite number of stationary solutions, the 
coarse-grained distribution f{r,\,t) can be trapped in one of them /q55(i", v) and remain 
frozen in that quasi stationary state until collisional effects finally come into play (on longer 
timescales). This steady solution is not always the most mixed state (it can be only partially 
mixed) so it may differ from Lynden-Bell's statistical prediction. Thus, for dynamical reasons, 
the system does not always explore the whole phase space ergodically. In general, the statistical 
theory of Lynden-Bell gives a relatively good first order prediction of the QSS without fitting 
parameter and is able to explain out-of-equilibrium phase transitions between different types 
of structures, depending on the values of the control parameters fixed by the initial condition. 
However, there are cases where the prediction does not work well (it can sometimes be very 
bad) because of incomplete relaxation. The difficulty is that we do not know a priori whether 
the prediction of Lynden-Bell will work or fail because this depends on the dynamics and it is 
difficult to know in advance if the system will mix well or not. Therefore, numerical simulations 
are necessary to determine how close to the Lynden-Bell distribution the system happens to 
be. Let us give some examples of complete and incomplete violent relaxation in stellar systems, 
2D turbulence and for the HMF model. 

3.2 Stellar systems 

The concept of violent relaxation was first introduced by Lynden-Bell [H] to explain the ap- 
parent regularity of elliptical galaxies in astrophysics. However, for 3D stellar systems the 
prediction of Lynden-Bell leads to density profiles whose mass is infinite (the density decreases 
as at large distances). In other words, there is no maximum entropy state at fixed mass 
and energy in an unbounded domain. Furthermore, it is known that the distribution functions 
(DP) of galaxies do not only depend on the energy e = + $(r) contrary to what is pre- 
dicted by the Lynden-Bell statistical theory. This means that other ingredients are necessary 
to understand their structure [16]. However, the approach of Lynden-Bell is able to explain 
why elliptical galaxies have an almost isothermal core. Indeed, it is able to justify a Boltzman- 
nian distribution / ~ e"^*^ in the core without recourse to collisions which operate on a much 
longer timescale t^eiax ~ {N/ In N)to [17\. By contrast, violent relaxation is incomplete in the 
halo. The concept of incomplete violent relaxation explains why galaxies are more confined 
than predicted by statistical mechanics (the density profile of elliptical galaxies decreases as 
instead of [18j). We note, however, that elliptical galaxies are not stellar polytropes so 
their DF cannot be fitted by the Tsallis distribution. 

For one dimensional self-gravitating systems, the Lynden-Bell entropy has a global maxi- 
mum at fixed mass and energy in an unbounded domain. Early simulations of the ID Vlasov- 
Poisson system starting from a water-bag initial condition have shown a relatively good agree- 
ment with the Lynden-Bell prediction [19j. In other cases, the Vlasov equation (and the corre- 
sponding A^-body system) can have a very complicated, non-ergodic, dynamics. For example, 
starting from an annulus in phase space, Mineau et al. [20] have observed the formation of 
phase-space holes which block the relaxation towards the Lynden-Bell distribution. In that 
case, the system does not even relax towards a stationary state of the Vlasov equation but 
develops everlasting oscillations. 

For three dimensional self-gravitating systems confined within a box, Taruya & Sakagami 
[21] found numerically that the transient stages of the collisional relaxation of the A^-stars sys- 
tem can be fitted by a sequence of polytropic (Tsallis) distributions with a time dependent q{t) 
index. Therefore, after the phase of violent relaxation, the system passes by a succession of 
quasi-stationary solutions of the Vlasov equation slowly evolving with time due to collisions (fi- 
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nite effects), until the gravothermal catastrophe associated with the Boltzmann distribution 
finally takes place. 

3.3 Two-dimensional vortices 

In the context of two-dimensional turbulence, Miller [22] and Robert & Sommeria [23] have 
developed a statistical mechanics of the 2D Euler equation which is similar to the Lynden- 
Bell theory (see [24j for a description of this analogy). This theory works relatively well to 
describe vortex merging [25] or the nonlinear development of the Kelvin-Helmholtz instability 
in a shear layer [SB]. It can account for the numerous bifurcations observed between different 
types of vortices (monopoles, dipoles, tripoles,...) [27] and is able to reproduce the structure of 
geophysical and jovian vortices like Jupiter's great red spot [281 129] . 

However, some cases of incomplete relaxation have been reported. For example, in the 
plasma experiment of Huang & DriscoU [30], the MRS statistical theory gives a reasonable 
prediction of the QSS without fit but the agreement is not perfect [31]. The observed central 
density is larger than predicted by theory and the tail decreases more rapidly than predicted 
by theory, i.e. the vortex is more confined. This is related to the fact that mixing is not 
very efficient in the core and in the tail of the distribution (these features can be explained by 
developing a kinetic theory of violent relaxation [32]). As observed by Boghosian ^33], the QSS 
can be fitted by a Tsalhs distribution where the density drops to zero at a finite distance. 

3.4 The HMF model 

The nature of the quasi stationary states (QSS) observed in the HMF model has generated an 
intense (and lively) debate in the community of statistical physics. 

Latora et al. [3lj performed A^-body numerical simulations starting from a water-bag initial 
condition with magnetization M(0) = 1 (unstationary). They observed the formation of QSS 
whose lifetime diverges with the system size A^. The Boltzmann distribution of statistical 
equilibrium fs = Ae~^^ is reached on a timescale treiax ~ N . In their Fig. 1(a), they compared 
the caloric curve T{U) of these QSS (bullets) with the caloric curve corresponding to the 
Boltzmann distribution (full line). They found a range of energies 0.5<t/<?7c = 3/4 where 
the caloric curve of the QSS disagrees with the caloric curve corresponding to the Boltzmann 
statistical equilibrium. We can interprete their results in another (complementary) manner. 
First of all, we note that the caloric curve of the QSS should be compared with the caloric 
curve predicted by the Lynden-Bell theory of violent relaxation since we are dealing with out- 
of-equilibrium structures (there is a priori no reason why the caloric curve of the QSS resulting 
from the violent collisionless relaxation should coincide with the caloric curve of the Boltzmann 
statistical equilibrium state resulting from the slow collisional relaxation). The question we now 
ask is: can the QSS be described by the Lynden-Bell theory? We note that, for a water-bag 
initial condition (two-levels), the distribution predicted by Lynden-Bell is similar to the Fermi- 
Dirac statistics f = 770/(1 + e^'^'^") [T5]. Furthermore, for the M(0) = 1 initial condition, we 
are in the dilute (non degenerate) limit of the Lynden-Bell statistical theory since the initial 
phase level r]o — >■ +00 [35]. Therefore, for that particular initial condition, we remark that the 
distribution predicted by Lynden-Bell coincides with the Boltzmann statistics f^^^ = Ae"^^, 
although it applies to the out-of-equilibrium QSS. Because of this coincidence, we can use 
the caloric curve T{U) reported in Fig. 1(a) of [33] to determine the domain of validity of 
the Lynden-Bell prediction. This curve shows that the Lynden-Bell prediction works well for 
U > Uc = 3/4 (i.e. in the region where the Lynden-Bell maximum entropy state is spatially 
homogeneous) and for U < 0.5 (i.e. in the region where the Lynden-Bell maximum entropy 
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state is strongly spatially inhomogeneous) . However, for 0.5 < U < Uc (i.e. close to the 
transition energy), the Lynden-Bell prediction fails. In that case, Latora et al. [31] show 
that the distribution has the tendency to remain spatially homogeneous {Mqss — 0) and that 
the velocity distribution is non-gaussian. These results strongly differ from the Lynden-Bell 
theory predicting a spatially inhomogeneous state with gaussian distribution /^^ = y4e~^^, the 
same as the statistical equilibrium state fs = Ae~^^. Therefore, close to the critical energy 
Uc, violent relaxation is incomplete and leads to a non-ergodic behaviour. Since standard 
statistical mechanics breaks down (standard statistical mechanics in the present context refers, 
in our sense, to the Lynden-Bell theory), Latora et al. [M] propose to describe this regime in 
terms of Tsallis generalized thermodynamics. This is an interesting idea to explore since there 
are no many other alternatives when the evolution is non-ergodic (another alternative could be 
to develop a kinetic theory of violent relaxation as attempted in [36]). 

Yamaguchi et al. [37] performed iV-body numerical simulations starting from a water-bag 
initial condition with magnetization M(0) = 0. For U = 0.69 > U* = 7/12, this initial 
condition is a stable steady state of the Vlasov equation. Furthermore, we remark that this 
spatially homogenenous water-bag initial condition is a minimum of energy E for a given mass 
M and phase level value 770 [35]. Therefore, this initial condition is the Lynden-Bell maximum 
entropy state. As a result, it does not evolve at all through the Vlasov equation. Yamaguchi 
et al. [37] show that it slowly evolves under the effect of collisions (finite effects) by passing 
through a series of stationary solutions of the Vlasov equation. This is similar to the results 
obtained by Taruya & Sakagami [21] for self-gravitating systems. Yamaguchi et al. [27] show 
that the Boltzmann distribution is reached on a timescale treiax ~ A^^'^^d and that the velocity 
distribution of the transient states is given by the curve reported on their Fig. 12. Recently, 
Campa et al. [38] obtained similar results and showed that these transient states can be fitted 
by a semi-elliptical distribution. Interestingly, we remark that a semi-elliptical distribution is 
a Tsallis distribution fg{v) = [fi — (3{q — l)t>^/2g]^/(''~^) with an index g = 3 (if we note 1 — g* 
instead of g — 1 this corresponds to g* = —1) □. Therefore, we observe that the results of 
Yamaguchi et al. [37] and Campa et al. [38], like the results of Taruya & Sakagami [21] in 
astrophysics, show that Tsallis distributions may he useful to describe the transient states of 
a collisional relaxation. In their paper, Yamaguchi et al. [37] (see also [39]) reject the Tsallis 
distributions because of the absence of power law tails in their curves of Fig. 12. However, 
power law tails are obtained only for a subclass of Tsallis distributions corresponding to indices 
g < 1 (in our notations [TOl SO]). For g > 1, the Tsallis distributions drop to zero at a finite 
value of the velocity, so they have a compact support. In particular, the distribution obtained by 
Yamaguchi et al. [37] and Campa et al. [38] appears to be well-fitted by a Tsallis distribution 
with g = 3 > 1 with a tail going to zero abruptly at a finite velocity Vmax- In view of the lively 
debate and the controversy about the applicability of the Tsallis statistics to the HMF model 
[STf [55] , it is amusing to realize that the distribution obtained numerically by Yamaguchi et al. 
[37] is in fact ... a Tsallis distribution! It has a compact support (g > 1) instead of power- law 
tails (g < 1). 

Antoniazzi et al. [41j performed A^-body numerical simulations starting from a water- 

^As a result, we can directly apply the stability criterion of [10] to obtain the critical energy above (below) 
which this distribution is Vlasov stable (unstable). Note first that the index q = 'd corresponds to a polytropic 
index n = 1 (see Eq. (144) of [TD]) or 7 = 2 (see Eq. (145) of [TO])- Therefore, according to the stability 
criterion (156) of [10], the critical energy is Ccrit = 1/7= 1/2- Using the relation U = e/4 + 1/2 (see Eq. (70) 
of [3S]) between the usual energy U used in [311 [371 [35] and the energy e used in [ini[3S], this leads to a critical 
energy Ucrit = 5/8. This coincides with the result obtained by Campa et al. |38j in a different manner. Note 
that the stability criterion (156) of [TUI which includes the Boltzmann (7 = 1), the Tsallis (7 = 1 + 1/n with 
n = 1/2 + l/(g — 1)), the water-bag (7 — 3) and the semi-elliptical (7 = 2) distributions is expressed very simply 
in terms of the polytropic index 7 similar to the one classically used in astrophysics [18] . 
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bag initial condition with energy U = 0.69 and magnetization M(0) between and 1. For 
this value of energy, the Lynden-Bell theory predicts an out-of-equilibrium phase transition 
from a homogeneous state to an inhomogeneous state above a critical magnetization Merit = 
0.897 discovered in [351 HI]. Numerical simulations [H] show that the Lynden-Bell prediction 
works relatively well for M(0) < Merit- This is confirmed by Campa et al. [2H] who find in 
addition that treiax ~ N^''^tr) as for M(0) = 0. However, above the critical magnetization, 
the results of Latora et al. [M] and Campa et al. [38] indicate that the Lynden-Bell theory 
does not work since the observed QSS is homogeneous {Mqss — 0) with non-gaussian tails 
while the Lynden-Bell theory predicts an inhomogeneous state {Mqss 0) with gaussian 
tails. This discrepency is particularly clear for the initial condition M(0) = 1 where the 
Lynden-Bell distribution coincides with the Boltzmann distribution f^^ ~ e"^*^ (non degenerate 
limit). Now, the early work of Latora et al. [31] indicates that this gaussian distribution is not 
observed and the recent work of Campa et al. [38] (for an isotropic water bag initial condition) 
shows that the QSS is well fitted by a semi-elliptical distribution. As we have seen, this is 
a particular Tsallis distribution with index g = 3 possessing a natural velocity cut-off. Such 
distributions, that rapidly drop to zero at a finite energy (here velocity) are typical products 
of incomplete relaxation. They are explained qualitatively by the fact that the high energy 
tail of the distribution in phase space does not mix well (a similar confinement is observed 
in the plasma experiment of Huang & DriscoU [30j discussed in Sec. 13. 3p . This confinement 
is consistent with a kinetic theory of incomplete violent relaxation [T6l [32| [36] . Therefore, as 
proposed in [25], the out-of-equilibrium phase transition predicted by the Lynden-Bell theory 
could be associated with a change of regime in the dynamics. For M(0) < M^rit, the system 
mixes well, the evolution is ergodic and the violent relaxation is complete leading to the spatially 
homogeneous Lynden-Bell distribution. Here, usual thermodynamics (in the sense of Lynden- 
Bell) applies. By contrast, for M(0) > M^rit, violent relaxation seems to be incomplete. In 
that case, the system does not mix sufficiently well, the evolution is non ergodic and the 
observed QSS differs from the Lynden-Bell prediction. This is associated with the appearance 
of fractal-like phase space structures, aging, glassy behaviour, power-law decay of correlations 
and anomalous diffusion. Rapisarda & Pluchino [42] have proposed to describe these features in 
terms of Tsallis thermodynamics. On the other hand, in a recent paper, Pluchino et al. [13] have 
shown explicitly that, in this non-ergodic regime, time averages and ensemble averages differ. 
The time averages can be fitted by g-distributions. We propose that the ensemble averages 
could also be fitted by a g-distribution with an index g > 1 leading to a natural velocity cut-off. 
As we have seen, the index g = 3 corresponds to the semi-elliptical distribution observed by 
Campa et al. [38]. Summarizing the above discussion, it seems that the Lynden-Bell prediction 
works relatively well far from the transition line separating homogeneous and inhomogeneous 
states in the Lynden-Bell theory (see the phase diagrams reported in [3S1 IH]) but that it fails 
close to this transition line: for fixed M(0) = 1 this is around Uc = 3/4 (Fig. 1(a) of [M] 
shows a discrepency with the Lynden-Bell theory in that region) and for fixed U = 0.69 this is 
around Merit = 0.897 (the simulations of [Ml 1121 |l3l |38] show a discrepency with the Lynden- 
Bell theory for M(0) ~ 1). In that case, we have non-gaussian velocity distributions, phase 
space structures, glassy dynamics, aging, anomalous diffusion 0... We must however be very 

^The kinetic theory developed by Bouchet & Dauxois [9 and Chavanis [TOl [4] is valid when the distribution 
of the bath is spatially homogeneous. When the velocity distribution has gaussian tails, like the Lynden-Bcll 
distributions obtained for < A/(0) < Alcrit 31] (the case M(0) = is special since the Lynden-Bell distribution 
coincides with the water-bag distribution with compact support), the velocity correlation function decays like 
{v{Q)v{t)) ~ {\nt)/t and the diffusion of angles is normal (with logarithmic corrections) [S]. If the velocity 
distribution of the bath is water-bag or semi-elliptic, like for A/(0) — [37l[38] or M(0) = 1 [Ml [38], standard 
kinetic theory [9l[4l[T0j predicts that the velocity correlation function has an exponential decay {v(0)v{t)) ~ e^*/"^ 



14 



careful because these striking features, like phase space structures and anomalous diffusion, 
could be due to finite size effects [16], HI] and disappear for +00 (note that N = 256000 
in [l6], N = 10^ in [41j while = 2000 in |42j). In the absence of phase space structures, 
we suggest that diffusion is normal because the bath distribution (semi-elliptical [38]) has a 
compact support (see footnote 5). 

Morita & Kaneko [47] performed A^-body numerical simulations starting from an initial 
condition with energy U = 0.69 and magnetization M(0) = 1 which is different from the water- 
bag. In that case, they find that the system does not relax to a QSS but exhibits oscillations 
whose duration diverges with (they find that the system relaxes towards the Boltzmann 
distribution on a timescale treiax ~ N). Therefore, the Lynden-Bell prediction clearly fails. 
This long-lasting periodic or quasi periodic collective motion appears through Hopf bifurcation 
and is due to the presence of clumps (high density regions) in phase space. We remark that 
this behaviour is relatively similar to the one reported by Mineau et al. [20] for self-gravitating 
systems, except that they observe phase space holes instead of phase space clumps. 



4 Slow collisional relaxation 

In this section, we discuss the process of slow collisional relaxation in relation with the kinetic 
equation (l33l) . 



4.1 About the i7-theorem 

When the system is spatially inhomogeneous, its collisional evolution can be very complicated 
and very little is known concerning kinetic equations of the form fl33l) . For example, it is not 
straightforward to prove by a direct calculation that Eq. fl5^ conserves the energy. However, 
since Eq. fl55]) is exact at order 0(1/A^), the energy must be conserved. Indeed, the integral 
constraints of the Hamiltonian system must be conserved at any order of the 1/A^ expansion 
(note that the neglect of collective effects in Eq. (l33l) may slightly alter the strict conservation 
of energy). On the other hand, we cannot establish the if-theorem for an equation of the form 
(I33p . It is only when additional approximations are implemented (markovian approximation 
and spatial homogeneity) that the if-theorem is obtained. To be more precise, let us compute 
the rate of change of the Boltzmann entropy 5*^ = — / ^ In ^dridvi with respect to the 
general kinetic equation fl33l) . After straightforward manipulations obtained by interchanging 
the indices 1 and 2, it can be put in the form 
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We note that its sign is not necessarily positive. This depends on the importance of memory 
effects. In general, the Markovian approximation is justified for A^ +00 because the correla- 
tions decay on a timescale Tcorr that is much smaller than the relaxation time treiax ~ Nto on 



leading to strictly normal diffusion of angles (this will be checked in a future contribution). However, if the 
system exhibits phase space structures like for M(0) > Merit and relatively small values of A'' [311112], the 
approach of Bouchet & Dauxois [9], which assumes spatial homogeneity, is not valid anymore. It is precisely 
the presence of these phase space structures that induces anomalous diffusion as studied in [42j . Therefore, 
there should not be any controversy since these authors [9] and [42] consider different situations as advocated 

in mm- 
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which the distribution changes (as discussed in Sec. 12.4^ this approximation is not completely 
obvious for self-gravitating systems and for systems that are close to the critical point). In that 
case, the entropy increases monotonically (see Sec. 12.31 for homogeneous systems). However, 
even if the energy is conserved and the entropy increases monotonically, it is not completely 
clear whether the general kinetic equation (155]) will relax towards the mean field Maxwell- 
Boltzmann distribution (1-24) of statistical equilibrium. It could be trapped in a steady state 
that is not the state of maximal entropy because there is no resonance anymore to drive the 
relaxation (this is the case for one dimensional homogeneous systems; see the discussion in 
Sec. 14.31) . It could also undergo everlasting oscillations without reaching a steady state. The 
kinetic equation (133|) may have a rich variety of behaviors and its complete study is of great 
complexity. 



4.2 The case of stellar systems 

The case of stellar systems is special. These systems are spatially inhomogeneous but, in order 
to evaluate the coUisional current in Eq. fl55]) . we can make a local approximation [IB] and work 
as if the system were homogeneous. This is justified by the divergence of the gravitational force 
F(2 1) when two particles approach each other so that the fluctuations of the gravitational 
force are dominated by the contribution of the nearest neighbour r2 — * ri [48j. The local 
approximation amounts to replacing f{Y2,^2,t — t) by /(ri,V2,t — r) in Eq. (155]) . This 
approximation is justified by the fact that the diffusion coefficient diverges logarithmically 
when Y2 — * ri (see below). Using the same argument, we can replace T^{2 — >■ 1) by F^{2 1) 
and JF^(1 ^ 2) by ^'^(l — > 2) = —F^{2 — > 1). We shall also make a markovian approximation 
/(ri, vi, t — r) ~ /(ri, vi, t), /(ri, V2, t — r) ^ /(ri, V2, t) and extend the time integration to 
infinity. Then, Eq. (155]) becomes 

Making a linear trajectory approximation Vj(t — r) = Vj(t) and rj(t — r) = rj — VjT, we can 
perform the integrations on ri and r like in Appendix A of Paper II. This yields the Vlasov- 
Landau equation 

ot OTi avi 
2.mGHnA^J,., (^___j/(.„v,t)/(r„v„t), (41) 

where InA = J^°°dk/k is the Coulombian factor ^18j. It must be regularized at small and 
large scales by introducing appropriate cut-offs, writing InA = \n{Lmax / Lmin) ■ Note that the 
divergence at large scales does not occur in Eq. (E^j. It only arises if we assume that the system 
is spatially homogeneous and infinite (and if we make a Markovian approximation and extend 
the time integral to infinity). In their stochastic approach, Chandrasekhar & von Neumann 
IH] argue that the Coulombian factor must be cut-off at the interparticle distance because the 
fluctuations of the gravitational force are described by the Holtzmark distribution (a particular 
Levy law) that is dominated by the contribution of the nearest neighbour. However, Cohen et 
al. [19] , considering a Coulombian plasma, argue that the integral must be cut-off at the Debye 
length, which is larger than the interparticle distance. This is confirmed by the kinetic theory 
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of Lenard [50] and Balescu [SI] which takes into account collective effects responsible for Debye 
shielding. In their kinetic theory, there is no divergence at large scales and the natural upper 
length scale appearing in the Coulombian factor is the Debye length. In plasmas, a charge is 
surrounded by a polarization cloud of opposite charges that diminishes the interaction. For 
gravitational systems, there is no shielding so we must stop the integration at R, the system 
size. Therefore, it is the finite spatial extent of the system that removes the Coulombian 
divergence. In a sense, the system size R (or the Jeans length) plays the role of the Debye 
length in plasma physics. On the other hand, the divergence at small scales comes from the 
break up of the linear trajectory approximation when two stars approach each other. This 
divergence also occurs in Eq. ( l33l) for the same reason: the unperturbed mean field motion 
( l32l) becomes incorrect when two stars approach each other. In fact, to obtain Eq. ( l33i) . we 
have assumed that the correlation function g is small with respect to /. This is true on average, 
but it is clear that correlations are important at small scales since two stars have the tendency 
to form a binary. Therefore, the expansion in powers ofl/N is not valid at any scale. One way 
to circumvent these difficulties is to use Eq. (155]) or (HT]) without modification but introduce a 
cut-off at the Landau length corresponding to a defiection of 90° of the particles' trajectory [§. 
Thus, we shall take ~ Clm/vlyp where Vtyp is the typical velocity of a star. Therefore, the 
Coulombian factor is estimated by In A = hi[Rv'^yp/Gm). Now, using a Virial type argument 
vlyp ~ (f^) ~ GM/R, we find that InA ~ InN. The relaxation time t^ due to encounters can 
be estimated from the Vlasov-Landau equation (j^Tj) by comparing the scaling of the l.h.s. and 



r.h.s. This yields l/t/j ~ mG'^lnAp/vf . The dynamical time is t^ ~ R/vtyp ~ l/\^pG where 



p ~ M/ R^ is the density. Comparing these two expressions, we get the scahng 

iH ~ ^to. (42) 

A more precise estimate of the relaxation time is given in [181 E] • The Vlasov-Landau equation 
conserves the mass, the energy (kinetic + potential) and monotonically increases the Boltz- 
mann entropy. The mean field Maxwell-Boltzmann distribution (1-24) is the only stationary 
solution of this equation (cancelling both the advective term and the collision term individ- 
ually). Therefore, the system tends to reach this distribution on a timescale {N/ In N)tD- 
However, there are two reasons why it cannot attain it: (i) Evaporation: when coupled to the 
gravitational Poisson equation, the mean field Maxwell-Boltzmann distribution (1-24) yields 
a density profile with infinite mass so there is no physical distribution of the form (1-24) in 
an infinite domain [Mt [55] . The system can increase the Boltzmann entropy indefinitely by 
evaporating. Therefore, the Vlasov-Landau equation fHTj) has no steady state with finite mass 
and the density profile tends to spread indefinitely, (ii) Gravothermal catastrophe: if the energy 
of the system is lower than the Antonov threshold Ec = —0.335GM^/R (where R is the system 
size), it will undergo core collapse. This is called gravothermal catastrophe because the system 
can increase the Boltzmann entropy indefinitely by contracting and overheating. This process 
usually dominates over evaporation and leads to the formation of binary stars [TH [55] . 



4.3 One dimensional systems 

One dimensional systems are also special. We have seen in Sec. 12.31 that one dimensional 
systems that are spatially homogeneous do not evolve at all on a timescale ~ NtD or larger 
because of the absence of resonances. However, if the system is spatially inhomogeneous, new 

^Note that the divergence at smaU scales does not occur in the binary encounter treatment of Chandrasekhar 
[52] and Rosenbluth et al. [53] which takes into account the exact two-body orbit of the particles instead of 
making a straight line approximation. 
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resonances can appear as described in [SS] so that an evolution is possible on a timescale A^^d- 
Then, we can expect that one dimensional inhomogeneous systems will tend to approach the 
Boltzmann distribution on the timescale Nto- To be more precise, let us consider the orbit- 
averaged- Fokker-Planck equation derived in [56]. Exploiting the timescale separation between 
the dynamical time and the relaxation time, we can average Eq. f l33l) over the orbits, assuming 
that at any stage of its evolution the system reaches a mechanical equilibrium on a short 
dynamical time. Therefore, the distribution function is a stationary solution of the Vlasov 
equation / ~ /(e, t) [where e = f ^/2 + $ is the individual energy] slowly evolving in time under 
the effect of "collisions" (= correlations due to finite N effects). Introducing angle-action 
variables, we get an equation of the form [56] : 

% = lijYl jrnA^n^'{lJ'?S{mn{J)-m'n{J'))[f{J')m^-f{J)m'^^ (43) 

m,m' ^ 

The important point to notice is that the evolution of the system is due to a condition of 
resonance between the pulsations fl{J) of the particles' orbits (this property probably extends 
to d dimensions but is technically more complicated to show). Only particles whose pulsations 
satisfy mfl{J) = m'Q{J') with (m, J) 7^ (m', J') participate to the diffusion current. This is 
similar to the collisional relaxation of two dimensional point vortices [57t [32] . It can be shown 
that Eq. fH3l) conserves mass and energy and monotonically increases entropy so that the 
system tends to approach the Boltzmann distribution of statistical equilibrium on a timescale 
~ Ntjj [Sn]- However, it may happen that there is not enough resonances so that the system 
can be trapped in a quasi stationary state different from the Boltzmann distribution. This 
happens when the condition of resonance cannot be satisfied so that mQ{J) 7^ m'Q{J') for all 
(m, J) 7^ (m', J'). In that case, the system is in a steady state of Eq. (H3l) which is not the 
Boltzmann distribution. This is what happens to point vortices in 2D hydrodynamics when the 
profile of angular velocity becomes monotonic [57]. In that case, the relaxation stops and the 
system will relax on a timescale larger than Nto- We may wonder whether the same situation 
can happen to systems described by a kinetic equation of the form (|33|) . 

5 Conclusions and perspectives 

In this paper, we have developed a kinetic theory for Hamiltonian systems with weak long- 
range interactions. A specificity of these systems is that they can be spatially inhomogeneous, 
which considerably complicates the kinetic theory. We have shown that the developement of 
correlations between particles creates a current in the r.h.s. of Eq. ([6]) that is the counterpart 
of the collision term in the Boltzmann equation for neutral gases. Therefore, for Hamiltonian 
systems with weak long-range interactions, the evolution beyond the Vlasov regime is driven by 
"correlations" due to finite effects. We have obtained a kinetic equation ([2SD valid at order 
0{1/N) that describes the evolution of the system on a timescale ~ Nti). For homogeneous 
systems, this equation reduces to the Landau equation. For d > 1, the Landau equation 
relaxes towards the Boltzmann distribution. Therefore, for > 1, the relaxation time scales 
like treiax ~ Nt^. This scaling has been predicted and observed in a spatially homogeneous 
two-dimensional Coulombian plasma [U EHl E] • This scaling probably remains true for spatially 
inhomogeneous systems in c? > 1 with the exception of self-gravitating systems that relax 
towards the mean-field Boltzmann distribution on a timescale treiax ~ {N/ In N)tD, unless they 
experience evaporation or gravothermal catastrophe. For one dimensional systems, like the 
HMF model, the situation is more complicated. For Vlasov-stable homogeneous systems, the 
kinetic equation (1331) reduces to df /dt = 0. Therefore, there is no evolution on a timescale of 
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the order ~ Nt^). We conclude that the relaxation time is larger than Nto- We could imagine 
that the evolution is due to three-body, four-body,... correlations leading to a relaxation time 
of the order of N'^t^, N^t £>■,■■■ However, Campa et al. [38], considering initial conditions with 
supercritical energy [/ > [4 = 3/4 for which the system is always spatially homogeneous, 
found that the relaxation time is extremely long scaling like t^eiax ~ e^. This suggests that the 
expansion of the BBGKY hierarchy in powers of may not be convergent in the homogeneous 
case and that another approach should be developed in that case. On the other hand, Morita & 
Kaneko [17] considering an initial condition with U < Uc and M(0) = 1, found a relaxation time 
of the order treiax ~ Nto- In their simulations, the system is always spatially inhomogenous 
(the magnetization in the oscillatory regime is non-zero). As explained in Sec. 14. 3[ spatial 
inhomogeneities can create new resonances that drive the relaxation towards the Boltzmann 
equilibrium (BE) on a timescale Ueiax ^ Ntc (predicted by the kinetic theory) that is much 
shorter than when the system remains spatially homogeneous. This could be an explanation 
(but not the only one) for the observed timescale in ^47] . Yamaguchi et al. [3^, considering a 
water-bag initial condition with U <Uc and M(0) = 0, found a relaxation time tj-eiax ~ N^tr) 
with 5 = 1.7. In their simulations, the system is spatially homogeneous but it progressively 
becomes Vlasov unstable and undergoes a dynamical phase transition from the homogeneous 
QSS to the inhomogeneous BE. This instability considerably accelerates the relaxation towards 
the Boltzmann equilibrium with respect to the case U > Uc ^38] where the homogeneous 
distribution remains Vlasov stable until the end (see below). The same phase transition happens 
in the simulations of Latora et al. [34J who considered a water-bag initial condition with U < Uc 
and M(0) = 1. Their system is roughly spatially homogeneous {Mqss — 0) but it also presents 
some phase space structures which may explain why they find a relaxation time of the order 
treiax ~ Nt^ shortcr than trdax ~ N^''^tD (we have seen that spatial inhomogeneities can 
accelerate the relaxation by creating new resonances). More generally, for water-bag initial 
conditions, it would be interesting to determine how the exponent 6 depends on the initial 
magnetization M(0) and energy U. Considering the phase diagram in (M(0), U) plane reported 
in [44J, we suggest that treiax ~ above the critical line Uc = 3/4 (no dynamical phase 
transition) and treiax ~ N^to below the critical line Uc = 3/4 (dynamical phase transition). 
For given initial magnetization M(0), we expect that 6 diverges as we approach the critical 
energy Uc = 3/4 above which treiax ~ e^. Below the critical line Uc = 3/4 and above the 
transition line Ucrit{M{0)) (the curve in Fig. 1 of [44J), observations [37] show that 6 ~ 1.7. 
Below the transition line Ucrit{M{0)), the QSS is either spatially inhomogeneous (according to 
Lynden-Bell's prediction if relaxation is complete) or spatially homogeneous with phase space 
structures in case of incomplete relaxation [311 112] • In that case, observations [31] show that 
5 ~ 1, which is consistent with the kinetic theory for inhomogeneous systems. Thus, for a 
given energy U < Uc, we expect that 6 passes from 5> 1.7 to 5 = 1 when M(0) overcomes 
the critical magnetization Mcrit{U) [33 IB]- These ideas, that are consistent with the partial 
numerical information that we have at present, demand to be developed in more detail. In 
fact, this is a first attempt to connect the relaxation time to the phase diagram [U, M(0)) and 
the situation may be more complicated than that. In particular, the relaxation time seems 
to strongly depend on the detailed structure of the inital condition. For example, using an 
isotropic water bag initial condition with M(0) = 1, Campa et al. [3H] find a relaxation time 
Ueiax ~ N^''^tD instead of the scaling t^ei ax ~ Nt£) reported by Latora et al. |34| . This may be 
related to the lack of phase space structures in the simulations of [38] . 

Finally, we conclude by proposing the following scenario similar to the one proposed by 
Taruya & Sakagami [2T] for self-gravitating systems. Analyzing the numerical results of [511 EZl 
[38] , we argue that, for many initial conditions, the transient states of the collisional relaxation 
of the HMF model can be described by spatially homogeneous Tsallis distributions (polytropes) 
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with a time varying index q{t) > 1 (i.e. n{t) > 1/2, 1 < 7(t) < 3) corresponding to a compact 
support. More precisely, we parametrize these transient states by a distribution function of 
the form f{v,t) = /(0,t)[l — f ^/fmaa;(^)^]"*-*''~^''^- It is easy to show that n{t) and Vmax{t) are 
related to each other by n{t) = 2TrVmax(ty / (kNe) — 1 where e = 4(f/ — 1/2) is the conserved 
energy and 2iT/{kN) = 1 in usual notations. This simple relation follows from Eqs. (28) and 
(29) of [10] for homogeneous systems where $ = 0, p = N/2it and p = E/n. It allows to 
determine n{t) by simply measuring Vmax{t)- These Tsallis distributions are quasi- stationary 
solutions of the Vlasov equation slowly evolving with time under the effect of collisions (finite 
N effects). Initially, n{t) is close to ra = 1 (i.e. 7 = 2, g = 3) corresponding to the semi- 
elliptical distribution observed by Campa et al. [38] in many circumstances. For U = 0.69 
(i.e. e = 0.76) and n = 1 we get Vmax — 1-23. Progressively, n{t) increases so as to attain 
the value n +00 (i.e. 7 = 1, g = 1) corresponding to the Boltzmann equilibrium state for 
t +00. For e < 1 (i.e. U < Uc = 3 / 4) , there exists a time t^, at which n(t) = ricrit = e/(l — e) 
(corresponding to 7(t) = 'jcrit = 1/e, lit) = qcrit = (e + l)/(3e — 1)) so that the Tsalhs 
distribution becomes Vlasov unstable (see the criterion (156) of [101) and the system rapidly 
relaxes towards the inhomogeneous Boltzmann distribution []. This accounts for the sudden 
dynamical phase transition observed in [27j from the homogeneous QSS {Mqss — 0) to the 
inhomogeneous BE (M = Mgg 7^ 0). For U = 0.69, the transition corresponds to ricru — 3.166 
leading to Vmax — 1-78 in qualitative agreement with [38]. In the supercritical case e > 1 
(i.e. U > Uc = 3/4), the Tsalhs distributions with q{t) > 1 are always Vlasov stable (since 
qcr^t = {e + l)/(3e - 1) = (4f/ - l)/(12f/ - 7) < 1 or, alternatively, e,„t = {q + l)/(3g - 1) < 1 
or Ucrit = (7g — l)/[4(3g — 1)] < Uc = 3/4) which explains the long lifetime behaviour observed 
by Campa et al. fSBI - This scenario suggests that Tsallis distributions can be attractors (or at 
least provide a good fit) for the transient states of the coUisional relaxation, for a large class 
of initial conditions {U = 0.69 < with M(0) = [371 [38] ; t/ = 0.69 < with M(0) = 1 
[5^ [S5] : and U > [3B])- Note, however, that the previous scenario is not valid for all initial 
conditions so these attractors are not universal. For example, in the numerical simulations of 
Antoniazzi et al. [H], the system relaxes towards a Lynden-Bell distribution with gaussian 
tails for < M < Merit = 0.897 and in the numerical simulations of Morita & Kaneko [^, the 
system develops everlasting oscillations. However, this picture now suggests to look in detail 
into the chaotic dynamics of the system in order to determine the basin of attraction of the 
Tsallis distributions [Ml [371 [2B| , the Lynden-Bell distributions [H] and the oscillatory states 
[U]. We hope to develop these issues, and check the above scenario, in future communications. 



A Some other kinetic equations 

If we make the Markovian approximation /(vi,t — r) ~ /(vi,t) and /(v2,t — r) ~ /(v2,t) in 
Eq. fl26l) but do not extend the time integral to infinity, we obtain 

^ = A^[ dY2K^''(w,t) /(Vi,t)/(v2,t), 

(44) 

with 

K^^iw, t) = (27r)'^m dr J dWk^'u^kf cos(k ■ wr). (45) 

^Note that Taruya & Sakagami [21 interprete this transition as a generahzed thermodynamical instabihty 
(in Tsalhs sense) while we interprete it as a dynamical instability with respect to the Vlasov equation J^0| /. 
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In (i = 3, the components of this tensor can be calculated by introducing a spherical system of 
coordinates with the 2;-axis in the direction of w. We find that 

K''''{w,t) = A{w,t) + B{w,t) —, 46 

yjA yjZ 



with 



, Svr^m /■+°° ,3 „ .,,,2 /■''"'V4sinr 4cosr\ , 
A{w,t) = / k^dku{k)^ ^ ] dr, (47) 



w Jo ' Jo \ 



„^ ^ IGn^m /'+°° ,3 „ .,,,2 f 2 sin T 4cosr 4sinr\ , 
B{w,t) = / k^dku{kf / + — dr. (48) 
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For t —>■ +00, the functions A and B reduce to 



5™ /"+00 
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A{w) = / k'uikydk, B{w) = 0, (49) 



and we recover the results (11-42) and (11-43) of Paper II. In o? = 2, the components of the 
tensor (H5i) can be calculated by introducing a polar system of coordinates with the x-axis in 
the direction of w. This leads to Eq. (H6|) with now 

A{w, t) = ^ edku{kf ^dT, (50) 

W ./n Jo r 



r+oo pkwt r 

'dr. (51) 



B(w,t) = ^^ / Pdkuikf / 
w Jo Jo 
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For t —>■ +00, the functions A and B reduce to 

A{w) = / k\{kfdk, B{w) = 0, (52) 

Jo 

and we recover the results (11-42) and (11-43) of Paper II. Finally, in c? = 1, we obtain 

K{w, t) = Anm e^kf^^^^^p^dk. (53) 

./n kw 



For t — i> -|-oo, we find that 

K{w) = 47i''mS{w) I ku{kfdk, (54) 



"+00 



which returns Eq. 0231) . Finally, we recall that Eq. fl26p ignores collective effects. As a simple 
generalization, we could replace in Eq. fl26|) the potential u{k) by the "screened" potential 
'u(A;)/|e(k, k- V2) | including the dielectric function (for the HMF model, this amounts to dividing 
the integrand of Eq. fl27|) by |e(l,i;2)P)- This leads to 

^ = (27r)^mA f dr [ d^^dm^ . ^,^ cos(k-wr) 
dt ov^ Jq J |e(k,k-V2)|^ 

xf7^-7^V(vi,t-r)/(v2,t-r). (55) 



dv'( 

Strictly speaking, this procedure is not rigorously justified for non Markovian systems since 
the dielectric function is obtained by assuming precisely that the distribution function does not 
change on the timescale of interest. Yet, this generalization could be performed heuristically in 
order to obtain a non Markovian kinetic equation taking into account some collective effects. 
If we make the Markovian approximation and extend the time integral to infinity, Eq. (l55l) 
returns the Lenard-Balescu equation. 
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